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ABSTRACT 



Breather solitons in a one-dimensional lattice of coupled nonlinear oscillators are 
numerically investigated. These are localized nonpropagating steady states that exist at 
frequencies either below the linear cutoff frequency (corresponding to the extended mode in 
which all the oscillators are in-phase) or above the upper linear cutoff frequency (corresponding 
to the extended mode in which each oscillator is 180° out-of-phase with its immediate neighbors). 
The lattice is damped and parametrically driven. A nonlinear Schrodinger theory, which assumes 
a modulational amplitude that is weakly nonlinear and slowly varying in space, is compared to 
the numerical data. The error is roughly 5 % at low amplitudes and 20% at high amplitudes. 
The regions in the drive parameter plane (amplitude vs frequency) where the breathers exist are 
numerically determined and compared to theory. A substantial discrepancy occurs at lower drive 
amplitudes where the theory predicts that the lower cutoff breather should exist, but where an 
instability is observed. Also in contrast to the theory, the region of the upper cutoff breather has 
relatively large areas in which quasiperiodicity occurs or the motion decays to rest. 
Quasiperiodicity is also observed in the lower cutoff breather. Finally, instead of a global 
parametric drive, an end drive is investigated. It is found that, for drive frequencies outside the 
linear propagation band, there is an amplitude threshold for the periodic ejection or "shedding" 
of propagating breather solitons. The quasiperiodicity that occurs for a global parametric drive 
may be a consequence of soliton shedding. 
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I. INTRODUCTION 



A. BACKGROUND 

A soliton is a localized wave whose shape or envelope is constant as a result of 
nonlinearities and dispersion. Furthermore, it collides elastically with other waves. Solitons 
were first observed as canal surface waves by John Scott Russell in 1834 (Dodd et al. 1978). 
Other systems are now known to support solitons of various types. Research in fiber optic 
solitons is currently very active because of their potential use in communications and optical 
switching (Mollenauer et al. 1991). 

In this thesis, a one-dimensional lattice of coupled nonlinear oscillators is numerically 
investigated. This system is known to possess kink solitons (Galvin 1990, Denardo et al. 1991). 
It will be shown that breather solitons can also exist. These are nonpropagating steady states in 
which most of the oscillators are approximately at rest while those in a localized region have 
substantial amplitude. Such a state is clearly not a linear mode of the system. Breathers are self- 
trapped states that occur at frequencies outside the linear propagation band. Two types exist: 
lower cutoff breathers, in which the oscillators are in-phase, and upper cutoff breathers, in which 
the oscillators are 180® out-of-phase. Lower and upper cutoff solitons can be considered as 
finite-amplitude modulations of the uniform lower and upper cutoff modes, respectively. 

In the absence of drive and dissipation, it is known that weakly nonlinear breathers 
described by the <#>“ theory are unstable, although the decay rate is extremely slow (Segur and 
Kruskal 1987). It is also known that instabilities can occur as a result of discreteness in systems 
that are undriven and undamped (Goedde 1990). In this thesis, dissipation and global drive are 
employed, which lead to stable breathers. 
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The simplicity of the model equation suggests that breathers can occur in a variety of 
systems. Indeed, the first observation of breathers was as cross surface waves on a long chatmel 
of liquid (Wu et al. 1984). Lower cutoff breathers have been observed in a pendulum lattice 
(Denardo 1990). Recently, upper cutoff breathers have been observed in a magnetically coupled 
pendulum lattice and confirmed numerically with an equation that models this system (Atchley 
1991). 

It is also shown that, for a local pure-ffequency drive confined to one end site of a lattice, 
it is possible to shed propagating solitons. The shedding occurs at a frequency that is 
quasiperiodic relative to the drive frequency. This phenomenon has previously been observed 
in only one system, surface waves just below the second cutoff mode in a large tank (Kit et al. 
1987, Shemer 1990). The existence of the soliton shedding in a simple lattice suggests that it is 
a general phenomenon which can occur in a variety of systems, and that a simple fundamental 
explanation should exist, although no such explanation is known at present. 

B. EQUATION OF MOTION 

To motivate the equation of motion that will be investigated, a model system is considered. 
The system is a uniform lattice of linearly coupled simple pendulums that oscillate transverse to 
the lattice. The torque on the nth pendulum is: 

M -1 ) " ^9^ sin(0„), J.B.l 

where fi is the torsional constant that characterizes the coupling, m is the pendulum mass, and 
1 is the pendulum length. For weakly nonlinear motion, the approximation 

sin(0) = e - - I.B.2 

6 
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I.B.2 



sin(0) = e - 1 ^ 

6 

can be made. Linear damping is included by adding to I.B. 1 a torque proportional to the velocity 
0'^ (with a negative coefficient), where the prime denotes differentiation with respect to time. 
To sustain localized structures in the lattice, a global parametric drive that results from vertically 
oscillating the lattice, is utilized. By the Equivalence Principle, in the reference frame of the 
lattice the effective acceleration due to gravity is: 

9eif “ ST + a cos(2o)t) , I.B. 3 



where a is the acceleration amplitude of the drive and 2o) is the drive frequency. 

Combining the above effects, and using Newton’s Second Law, the equation of motion of 
for the lattice can be expressed as: 

- 20„ + 0„_i) +p0'n+ [o)„2+2ncos (2(Ot) ] 0^ - a0^^, I.B. 4 

where c^ is a measure of the coupling between lattice sites and w, is the cutoff frequency of a 
pendulum. For simplicity, the term proportional to the product of the drive and nonlinearity has 
been dropped. It can be shown that this has no essential effect upon the localized structures. The 
nonlinear coefficient a equals a)„^/6 if (1.B.4) is to approximate a pendulum lattice. The model 
is generalized by allowing the nonlinear coefficient to be positive or negative. The lattice is 
driven parametrically with an amplitude of 2t). The fundamental response frequency w is half the 
drive frequency. The lower cutoff mode of the lattice is characterized by uniform motion in the 
lattice (Fig. I.B.l), and has linear frequency The upper cutoff mode has a 180° phase 
difference between adjacent lattice sites (Fig. LB. 2), and has linear frequency (oj,^ + 4c^*'^. For 
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a softening (a > 0) lattice, the upper cutoff mode is stable. The lower cutoff mode is subject 
to the Benjamin-Feir instability (Denardo 1990), and the motion evolves into one or more 
breather solitons (Fig. I.B.3). For a hardening (a < 0) lattice, the roles of the cutoff modes are 
reversed and an upper cutoff breather evolves (Fig. I.B.4). 
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II. THEORY 



A. LOWER CUTOFF BREATHER 

For modulations of the lower cutoff mode, the displacement (6J in I.B.4 can be written 
approximately as: 

d t) complex conjugate, ll.A.l 

where /l(x,t) is a complex differentiable function of its arguments. The lattice spacing is assumed 
to be unity. A displacement (0J containing a third harmonic term produced approximately the 
same results as ll.A.l and was consequently discarded. Substituting ll.A.l into I.B.4 and 
requiring that A be slowly varying and weakly nonlinear, one obtains the modulational equation 

2i(tiA - c^ A^ + - 0 )^ + i(oP)A + r)A* - 2a \ A^ A, ll.A.l 

which is a nonlinear Schrodinger equation. If the modulation is assumed to be nonpropagating, 
a stationary solution 

A (x,t) - A (x) II. A. 3 

exists where A and 5 are real. Substituting II. A. 3 into II. A. 2 gives: 

- A^' + i(A> - 3a + iti)PA - II. A. 4 

Equating the imaginary parts of both sides yields an expression for 6: 

sin(26) - . II. A. 5 

n 
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The real part of II. A. 4 yields: 



- \l A + 3 a - 0 , II. A. 6 

where \i - (ji - (j>^ + \jr\^ - . 



There are actually two possibilities because of the radical. However, for a softening lattice the 
negative branch leads to an unstable solution (Denardo 1990). 

Equation II. A. 6 can be integrated to yield an elliptic function. A localized solution to 
II. A. 6 is of the form A = a sech(bx), provided 



a - 



2c^ 
N 3a 



and b - 



N 



JL 

/n 2 



Then 



A - 



££■ sech 




{x - xJ 


3a 




J 



II. A. 1 



Substituting the solution for A into II. A. 1, the expression for the displacement is: 



e, 






JM, 

3a 



sech 



^ (x - xj cos((i)t + 6) . 

C2 



II. A. 8 



This solution is valid for /r > 0 and a > 0 (softening). 

The physical reasoning for the existence of the breather is derived from the curvature of 
the modulation envelope. For the softening lower cutoff breather, negative curvature in the body 
of the envelope produces a restoring force (Fig. II. A. 1). The frequency in the body is below 
cutoff due to the nonlinearity which decreases the frequency more than the curvature increases 
it. In the tail of the envelope the curvature is positive, thus an anti-restoring force. The 
frequency is below cutoff and the breather is evanescent in the tail. At the inflection points, the 
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nonlinearity is the sole reason the frequency is below cutoff. The breather is thus a self-trapped 
state due to the combination of the nonlinearity and curvature of the modulation. As a result, 
steady state motion is possible. This motion has been observed in the simulated lattice. 

Concluding that steady state motion is possible, one should be able to determine the values 
of Ti and CO where a stable solution exists. Requiring fi from II.A.6 to be positive, 

Tl^ > ( 0)^2 _ )2 + ^^2 p2^ ' II.A.9 

If a) = o)„, this yields the Matthieu hyperbola. The sech solution (Eqn. II.A.8) is unstable inside 
the hyperbola because of growth in the wings. There are two conditions for the solution to exist: 
j; > CO/S and n > 0. The first condition is ij > cojS if co»co„. The second condition is related 
to II.A.9, which can be rewritten as 

V Tl^ - 0)^ _ (^2|_ JJ.A.IO 

If CO > co„ then /c > 0. However, if co < co„ the requirement for /r > 0 is i; > co/3. The sech 
solution decays to rest below the boundary i] > co/S. 
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B. UPPER CUTOFF BREATHER 



The displacement (6J for modulations of the upper cutoff mode can be written as 

6^- (-1) ” A (u, t) complex conjugate II.B.l 

Where A(n,t) is a complex differentiable function of its arguments. Again, the higher order terms 
have been neglected, and the lattice spacing is assumed to be unity. The equation that describes 
the weakly nonlinear complex amplitude is again a nonlinear Schrodinger equation: 

2io)Aj. + - 0 )^ + io)P)A + T| A - 3a I A P A. II. B. 2 

Equation II. B. 2 is identical to II. A. 2 with the substitutions and c^-*-c^, where co,^ = 

+ 4c^, which is the square of the linear frequency of the upper cutoff mode. Assuming a 
nonpropagating modulation, a stationary solution for A is: 

A {x,t) - A (x) II. B. 2 

where A(x) and 5 are real. Substituting II. B. 3 into II. B. 2 gives: 

a'' + - 0 ) 2 )^ - 3a A^ + ioPA - -tiAe‘2^*. JJ.B.4 

Setting the imaginary parts of both sides equal and solving for 6, the phase relation is again: 

sin(26) - II. B. 5 

n 

The real part of II. B. 4 yields: 

-c^ A^^ + vA + 3aA-^-0 II. B. 6 

Where v - - w^p^ 
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which is similar to I1.A.6. Again, the selected sign of the radical leads to a stable solution. A 



trial solution of the form A = a sech(bx) used in II.B.4 yields: 






II.B.l 



Substituting the solution for A into II.B.l, the displacement (0J for the upper cutoff mode is: 



which is valid for »'>0 and a<0 (hardening). 

The physical mechanism for the existence of the upper cutoff breather is essentially 
identical to that of the lower cutoff breather. The negative curvature in the body of the hardening 
upper cutoff breather now produces an anti-restoring force (Fig. II.B.l) (Denardo 1990), 
however, the nonlinearity is stronger than the curvature and the frequency in the body is above 
the linear upper cutoff frequency. The positive curvature in the tails of the modulation creates 
a restoring force and the frequency is again above the cutoff frequency. Similar to before, the 
frequency is above the cutoff frequency at the inflection points as a sole result of the nonlinearity. 
Thus, the upper cutoff breather is also a self-trapped state due to the combination of curvature 
and nonlinearity. 

As for the lower cutoff case, a region of stable, steady state motion should also exist for 
the upper cutoff breather. Constraining p from II. B. 6 to be positive. 




II. B. 8 



r\^ > ( 0)^^ - 0)^ )^ + 0)^ p^. 



II. B. 9 
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Again, if w » a;„ this yields the Matthieu hyperbola. The wings of the sech solution are unstable 
inside the hyperbola. There are two conditions for the solution exist: 17 > co/3 and 1 / > 0. The 
first condition is j; > cojS if co*co„. The second condition is related to II. A. 20 which can be 
rewritten as 



\J r\^ - > IcOj^^ - <1)^1. JJ.S.IO 

If CO > coi, then n > 0. However, if co < coj the requirement for /t > 0 is 17 > co/3. The sech 
solution decays to rest below the boundary 17 > co/3. 
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C. END DRIVEN LATTICES 



Nonlinear lattices driven at one end and at a frequency outside the linear propagation band 
have been observed to periodically eject solitons (Sect. III. D). This is in strong contrast to the 
linear theory which predicts only a steady state evanescent wave, as will now be shown. The 
equation of motion for a lattice of linearly coupled nonlinear oscillators (I.B.4) can be modified 
to replace the parametric drive with a direct drive. Further restricting the drive to the end site 
alone and removing the periodic boundary condition yields: 

73 - 0 : 0o-C^(9^ - 0Q ) +P0'o + <O 0Q + T1COS (q t) - a0Q^ II. C. 2 

n>o-. e"„-cue„„ - e„ - «e„=. 

This produces an end driven lattice. The lattice is still linearly coupled and nonlinear, but is now 
more suitable for the observation of propagating solitons which have a spatially varying phase 
(Denardo 1990). Focusing on the linear motion in the continuum limit, equations II. C. 2 have 
the form; 

X > 0 : - c2 + p + 0)^2 y - 0 II. C. 3 

o / \ 

X - 0 : Vx ^ o y ^ ^ cos (cot) . 

3 . 

Consider a solution of the form: 

y - A + complex conjugate, II. C. 4 

with a complex wavenumber x = k + ia. Substituting II.C.4 into II. C. 3 yields: 

-0)2 - (iic - a)^ - iwP + o)^^ - 0. II. C. 5 



14 



Separating the real and imaginary parts of II. C. 5 gives: 



Re: 



a 



2 _ . 



(0 



- 



Im: 



a k 



2c2 ■ 



The case of interest is w < Wq. The solution for )S=0 is: 



k-0 



and 




II. C. 6 



Physically, the wave is purely exponential in space; the sign of the displacement is the same for 
every lattice site. If 0 then k can be eliminated using the real and imaginary parts of II.C.5 
to obtain: 



a* - 2 



- 0)2 

- 



2c‘ 



G)P 

2c2 



- 0 . 



II. C.l 



Solving for and choosing the positive root , since a was assumed to be real, gives: 



[o)„2 - (j2 + ^(<0^2 _ ^2)2 ^ 2j_ II.C.B 

The solution for k is then: 

Ic . A . II, C.9 

2c 2 a 

If the solutions for oP and k are specialized to weak damping, w/S « then: 



a - 






- 



0)' 



and k 



ojP 

2c 



v/oT 



II. C. 10 



- U)‘ 
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The damping gives rise to a spatial phase (k 0). Specializing to strong damping, w/S > 
or near-propagation, yields: 



a 




JJ.C. 11 



This is the largest that k/a can be (unity). In general, O^k/a^ 1 (Fig. II.C.l). 
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Figure II.C.l Evanescent Wave For k/a=l 



17 



III. NUMERICAL SIMULATION 



A. IMPLEMENTATION 

The motion of the lattice of linearly coupled nonlinear oscillators was solved with the 
Euler-Cromer method. This implementation utilized acceleration as the basic numerical iteration. 
The position calculation of each lattice site after a given time step was divided into three sections: 

1 . The calculation of acceleration for each lattice site due to the position of it and the 
two adjacent sites. 

2. The calculation of the new velocity for each site by multiplying the acceleration by the time 
step, and adding the old velocity. 

3. The calculation of the new position for each lattice site using the new velocity for each site. 
Periodic boundary conditions were imposed, thus making the lattice a ring lattice with any 
curvature ignored. The Euler-Cromer method converged, with an increasing number of time 
steps, to a "true" lattice amplitude (Fig. III. A. 1). A time step that gave an amplitude within . 1 % 
of the "true" response amplitude allowed the program to run at a fairly high rate of speed with 
minimal loss of accuracy. 

The program is highly interactive and gives a "real time" picture of the lattice motion. The 
numerical simulation imposes some interactive restrictions. The elapsed time in the model frame 
must be reset periodically (preferably every period of the drive) to avoid strong transients when 
changing the drive frequency. The problem stems from the parametric drive term (cos(2o)t)) 
since the time increment is based on the response period. Thus a change in frequency in the 
middle of the cycle produces a discontinuity in the size of the time increment. The solution is 
to employ an integer counter to reset the elapsed time after one period of the drive and to require 
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any frequency changes to occur when the time is reset. The startup of a modulation envelope or 
"profile" can be accomplished in two ways: restarting a saved profile, or starting a profile created 
from theory. A profile can be saved while the program is running, thus eliminating the need to 
duplicate steps required to reach a particular set of parameters. All necessary parameters such 
as drive amplitude and frequency in addition to the position and velocity for each lattice site are 
written to a file of the researcher’s choice. The modulation must meet certain stability criterion 
to be recorded and is "captured" at the upper turning point (response) of a selected site. A 
previously saved profile can be recalled using only the filename. Since the response and drive 
have a phase difference, the program restarts the profile with a time phase correction to minimize 
startup transients. A startup file could also be generated using one of the programs that calculate 
a profile from theory using keyboard input of drive, coupling and other parameters. The 
equations of motion are in a subroutine, allowing them to be altered or replaced quickly. A 
different method of calculating the lattice site positions would probably not cause any difficulty 
as long as the new coordinates were available at the same place in the program as before. A 
complete list and more detailed description of the program options as well as a program code 
listing are included in the Appendix. 
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B. LOWER CUTOFF BREATHERS 



In order to investigate the stability of lower cutoff breathers, a drive plane, consisting of 
drive amplitudes and frequencies where the particular structure exists, was obtained by numerical 
solution of the equation of motion. The structure chosen for the lower cutoff drive plane was a 
symmetric breather (Fig. IlI.B.l) in which the maximum amplitude is centered on a lattice site 
(on-site). The nonlinear parameter (a) was 1 (lower cutoff), the coupling (c^ or 7) was .1 
(weak), dissipation (fi) was .03 (small) and the lattice consisted of 50 sites long. The points in 
the drive plane (Fig. III.B.2) were almost exclusively obtained by incrementing the drive 
amplitude by .1% (or less) and recording the amplitude at which transitions occurred. 
Additionally, to probe stability, a <+3% (of the maximum response amplitude in the lattice) 
random kick was applied to the lattice after each increment. 

The upper boundary from w=.76 to .89 was characterized by the lattice "going over the 
top" after perturbation (random kick). Because, the potential well for the lower cutoff mode 
(Fig. III.B.3) is not infinitely high, the lattice is able to escape the well (go over the top) and 
become unstable. There were no visual indications in the lattice of an impending instability 
before the boundary was reached. The right upper boundary, from w=.89 to .94, was marked 
by the generation of another symmetric breather (out of phase with the first) on the opposite side 
of the lattice ring from the initial breather. The lattice then went over the top from a random site 
fairly quickly after the generated breather reached an amplitude comparable to the original 
breather. Further tests with a lattice of 100 sites indicated that the generated breather is 
approximately 27 sites away from the center of the original breather and will arise on both sides 
of the original if the sites are available. The breather structure in the extreme lower right corner 
is quasiperiodic until the boundaries where it goes over the top without the generation of another 
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breather. The lower boundary from w = .85 to .947 is preceded by a quasiperiodic region. 
Ultimately the lattice either goes over the top or decays to rest after very strong transients. 
Careful observations are suggestive of quasiperiodicity resulting from the energy loss due to 
soliton shedding. The left lower boundary from «=.76 to .85 was only preceded by the 
quasiperiodic region at the points indicated but the strong transients were consistently observed 
before the lattice went over the top. A comparison of the stable region predicted by theory for 
the lower cutoff symmetric breather shows good agreement for the upper right boundary but a 
discrepancy on the lower boundary (Fig. III.B.4). Some destabilizing factor has prevented the 
lower boundary from reaching the theoretical limit. Experimentation with a half-site centered 
symmetric breather (Fig. III.B.5) revealed a link between the strong transients observed in the 
on-site breather drive plane but did not extend the lower boundary to the theoretical limit. The 
strong transients occur at the same drive amplitudes that the half-site breather shows quasiperiodic 
behavior. Increasing the drive amplitude results in a transition from half-site to on-site at 
approximately the start of the quasiperiodic line in the on-site drive plane. The half-site breather 
is stable as drive amplitude is decreased, but it decays to rest at approximately j ]= 0.9 for «=0.9. 
The transition from steady state motion to rest is interesting in that it shifts back to an on-site 
breather before it decays to rest. 

Comparison of the theoretical profile for the parameters used in the drive plane was done 
for the lowest and highest amplitude saved profiles as well as a half-site profile. The lowest 
amplitude comparison was sharper and lower in amplitude than the predicted profile (Fig. 
III.B.6). The highest amplitude comparison was also sharper and even lower in comparable 
amplitude (Fig. III.B.7). A comparison of the half-site profile shows much closer agreement with 
theory although the data still exhibits a greater amplitude (Fig. II1.B.8). The amplitude of the 
saved profile from the program was not exactly the maximum amplitude of that particular state 
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did not precisely coincide with the peak response amplitude. This was corrected by using the 
maximum coordinate plus the coordinates for the preceding and following time steps to find the 
maximum amplitude by parabolic curve fitting. By recording the time step, the data could be 
scaled according to the time difference between the turning point and actual capture time. 



23 




3pnt).i I du^ 



Figure III.B.l Lower Cutoff Breather 
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Figure III.B.2 Lower Cutoff Drive Plane 
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Figure III.B.3 Lower Cutoff Potential Well 
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Figure III.B.4 Lower Cutoff Drive Plane and Theory 
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Figure III.B.5 Symmetric Half-Site Lower Cutoff Breather 
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Figure III.B.6 Lower Cutoff Breather Theory Comparison 
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Lower Cutofr Breather Comparison To Theory 




Theopij: line Data: «» Lattice Site eta=.099 omega=.77e 

Figure III.B.7 Lower Cutoff Breather Theory Comparison 



Lower Cutoff Breather Comparison To Theory 
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Figure III.B.8 Lower Cutoff Breather Theory Comparison 
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C. UPPER CUTOFF BREATHERS 



The structure chosen for the upper cutoff drive plane was the symmetric on-site breather 
(Fig. III.C.l). The parameters for the drive plane were a=-l (hardening), 7 =.l (weak 
coupling), /S=.03 (small) and a lattice length of 50 sites. The points in the drive plane (Fig. 
III.C.2) were almost exclusively obtained by incrementing the drive amplitude by .001 and 
recording the amplitude at which instabilities were observed. After each increment, a random 
kick of < +3% (of the maximum lattice response amplitude) was applied to the lattice. 

The upper left boundary marks the sharp transition from a breather structure to a complex 
upper cutoff mode structure. The transition time was fairly short and was characterized by the 
growth of one section of the lattice to an amplitude comparable to that of the original breather. 
The growing section was in a pure upper cutoff mode and spread to include the rest of the lattice 
after reaching full amplitude. The area inside the circles (Q) was where the lattice motion was 
quasiperiodic. The boundary was difficult to ascertain for omega 1.25 and higher since local 
regions appeared to exhibit quasiperiodicity but in fact merely had long settling times. Regions 
marked H denote areas of no lattice motion. Inside the region and below the lower boundary, 
the lattice decayed to rest with minimal transients. Region QS was characterized by a 
quasiperiodic shedding breather. The quasiperiodic breather was usually observed near the 
boundaries of QC and was similar to Fig. III.C.l. The center site of the breather was 
quasiperiodic and small solitons were ejected fi'om the center in both directions simultaneously. 
The quasiperiodic breather spontaneously transitioned to a complex cyclic state at the lower 
boundary of QC (asterisks). The lattice motion in the small asterisked area (QC) was of two 
general types, a complex cyclic state or an anti-symmetric half-site breather. The cyclic state that 
evolved from the quasiperiodic breather appeared to be chaotic and then settled into multiple on- 
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site breathers. The motion shifted back and forth between the two states with a fairly long cycle. 
If the lattice was in this particular mode when the boundary to H (right side of QC) was reached, 
the cycle was broken and the lattice settled into a multiple breather state. The multiple breathers 
merged two by two until a single on-site breather was left. The single breather then decayed to 
rest. The lattice also spontaneously shifted from the complex cycle to a half-site antisymmetric 
breather (Fig. III.C.3). A half-site breather consistently evolved from an quasiperiodic on-site 
breather at the upper boundary of QC. The half-site breather exhibited a degeneration similar 
to that of the complex cyclic at the boundary to region H. The breather spontaneously changed 
back into an on-site breather and then decayed to rest. Multiple antisymmetric half-site breathers 
were also observed. A state with two half-sites 10-15 sites apart (center-to-center) were believed 
to be transferring energy by soliton shedding. A growth in amplitude in the center sites of one 
breather corresponded to the arrival of a soliton from the other breather and the expected 
decrease in the center-site amplitude of the second breather due to soliton ejection was observed. 

The upper and lower boundaries compare well with the theoretical boundaries (Fig. 
III.C.4). The upper boundary does not quite follow the hyperbola at the higher drive amplitudes, 
however the theory is approximate at these amplitudes. The lower boundary, which corresponds 
to the lattice decaying to rest, coincides almost exactly with the predicted threshold of lattice 
motion. Comparison of theoretical modulation envelopes and the numerical data were 
completed for extremely high, high and low response amplitude cases as well as a half-site case. 
The observed amplitudes were adjusted to that of the turning point as described in Section II. B. 
The normal 180® phase difference between lattice sites was eliminated for easier modulation 
envelope comparison. The extreme case was not as sharp as predicted and was »70% of the 
theoretical amplitude (Fig. III.C.5). Similarly, the amplitude of the high case was *90% of 
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theoretical but the envelope was sharper than expected (Fig. III.C.6). The low amplitude case 
exhibited the inverse characteristics of the extreme case (Fig. III.C.7). The observed envelope 
was sharper than the theory and *=5% higher in amplitude. A comparison of the half-site case 
showed very good agreement with theory (Fig. III.C.8). 
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Figure III.C.l On-Site Upper Cutoff Breather 
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Figure III.C.2 Upper Cutoff Drive Plane 
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Figure III.C.3 Anti -Symmetric Half-Site Upper Cutoff Breather 
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Figure III.C.4 Upper Cutoff Drive Plane and Theory 
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Figure III.C.5 Upper Cutoff Breather Theory Comparison 



Upper Cutoff Breather Conparison To Theory 
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Figure III.C.6 Upper Cutoff Breather Theory Comparison 
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Upper Cutoff Breather ComparisoTi To Theory 
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Figure III.C.7 Upper Cutoff Breather Theory Comparison 
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Figure III.C.8 Upper Cutoff Breather Theory Comparison 
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D. SOLITON SHEDDING 



The initial states for both upper and lower cutoff mode response investigations were 
evanescent states in a lattice of 100 sites. One end of the lattice was driven and the other end 
was given a free boundary condition. The large number of sites minimized reflection while 
maintaining the speed of the simulation. The parameters for the upper cutoff were 7=. 75, 
/S=.03, while the lower cutoff were 7=. 50, /5=.03. Snapshots in time of the on-screen motion 
of the lattice were obtained to help describe the motion. In these "snapshot" figures, each dot 
is an individual lattice site and the site furthest to the left is the driven end site. 

Soliton shedding was observed in both the upper and lower cutoff, end driven lattices. The 
shedding was similar for both modes, with differences in the actual method of energy transfer 
from the driven site. The upper cutoff shedding was produced by the cyclic relaxation in 
amplitude of the driven end site. The lower cutoff shedding was also generated by end site 
relaxation but transferred energy through a "pivot" point before the soliton was created. 

An artificial potential well of the form: 

f(x) — 

V^l + 

was used to replace f(x)= -x + x^ in the lower cutoff program to permit amplitudes of the 
magnitude necessary to observe shedding in the softening mode. Otherwise, the oscillators would 
"go over the top." The upper threshold was marked by a distinct motion change and soliton 
shedding (Fig. III.D. 1). The shed solitons were very small and strongly evanescent from co = .95 
down to o) = .6. The soliton shedding at o) = .6 was distinct and the structure easily visible until 
20-30 sites from the driven end. The shedding results were similar for co= .5 with solitons visible 
until approximately site 35. A third harmonic component appeared at low drive amplitudes 
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(r;= . 144) and continued until the shedding threshold. The harmonic was identified by taking the 
FFT of a time series from the driven site (Fig. III.D.2). The visible result of the harmonic was 
an evanescent ripple wave (Fig. III.D.3). 

A fifth harmonic was identified for w=.3 (Fig. III.D.4) and it also disappeared at the 
shedding threshold. The strength of the harmonic appeared to vary slowly with amplitude but 
over a fairly large range. The driven end site became quasiperiodic from 7; = .415 to 1.42 and 
from 2. 14 to the shedding threshold for w= .2. Small wave packet ejection was observed from 
71=2.00 to 2.59. Upon reaching the shedding threshold, the end site relaxed and a large soliton 
was shed (Fig. III.D.5). The soliton would propagate 6-8 sites, then appear to shed a much 
smaller soliton. Both the original and shed soliton would then damp out quickly (Fig. III.D.6). 
A seventh harmonic appeared for «=.2 but was not visually identifiable. Odd harmonics up to 
the 11th were identified for o)=.l (Fig. III. D. 7-8) and visually affected the lattice motion. The 
end site quasiperiodicity started at r}=2.25 and the evanescent wave packets started at tj= 1.35. 
The shedding threshold was not reached by tj=6.0. 

The dynamic threshold for the upper cutoff lattice exhibits a consistent trend even though 
the motion of the resultant state changes significantly between tj= 2.32 and 2.33 (Fig. III.D.9). 
The basic motion changes from soliton shedding to a non-shedding evanescent state. Shedding 
for o)=2.025 occurs very slowly and the soliton propagates 30-40 lattice sites and then appears 
to stop (Fig. III.D. 10). Apparently the nonlinearity lowers the frequency enough that linear wave 
propagation can occur inside the structure and it appears to "ring" as it decays (Fig. III.D. 11). 
The shedding occurs slowly for a>=2.05 and the soliton is large but jumbled (Fig. III.D. 12) and 
dies out quickly. The size of the shed solitons is cyclic for a>=2. 1-2.3 and the structure is very 
diffuse (Fig. III.D. 13). Several small solitons are ejected followed by one large one and the 
relaxation of the end site is very distinct. Intermittent shedding occurs for w=2.31 and the first 
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8-9 sites move together (Fig. III. D. 14). The transition point from shedding to nonshedding 
seems to be between w = 2.32 and 2.33. The motion for co=2.32 is varied in the number of sites 
in tandem, the size of the shed solitons and the magnitude of the end site relaxation. The 
shedding is transitional for w=2.33 and the motion settles into a complex evanescent state. The 
first 6-9 sites move in tandem and the last site is the first of 2-3 evanescent sites. The motion 
for co = 2.35-2.6 is similar but only 2-3 sites are in tandem. The vertical boundary between 
shedding and non-shedding states was confirmed by obtaining high drive level evanescent states 
and crossing the boundary into the shedding region. Specifically, the states co=2.35 tj= 5.05 and 
«=2.5 7 j = 8.0 still exhibited the same lattice motion initiatated at the transition boundary. The 
states were then decremented in frequency to cross into the shedding region. The motion of both 
states slowly evolved to the previously described lattice motion for the final frequency. 
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Figure III.D.l End-Driven Lower Cutoff Drive Plane 
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Figure III.D.2 Lower Cutoff FFT Spectrum w=.5 




Spectrun for elenent O LATTICZX.C 

Alpha: 1.000000 Beta: 0.030000 Ganna : 0.500000 Eta: O.GOOOOO Onega: 0.300000 
rtax anp = 596.79-1171 Itax frcq shown: 1 . 19366Z Tine interval: 0.Z09-1-10 




Figure III.D.4 Lower Cutoff FFT Spectrum w=.3 
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Figure III.D.7 Lower Cutoff FFT Spectrum w=.l 
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Figure III.D.9 End-Driven Upper Cutoff Drive Plane 



46 




Eta= . 81 Omega=2 . 025 



*• 



Figure III. D. 11 Time Snapshot o)=2.025 
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IV. SUMMARY AND CONCLUSIONS 



It has been shown numerically that nonpropagating steady state breather solitons can occur 
in one-dimensional lattices of coupled nonlinear oscillators that are damped and parametrically 
driven. The simplicity of the model suggests that these states can exist in a variety of systems. 
Amplitude data agree well with a nonlinear Schrodinger theory at low amplitudes, but disagree 
substantially at higher amplitudes. The theory assumes that the amplitudes are weakly nonlinear 
and slowly varying in space. The continued existence of breathers at amplitudes where the theory 
breaks down is evidence of the robustness of these states. Dissipation and drive stabilize the 
breathers, although this is not yet understood. 

Breathers exist for a range of drive parameters (amplitude and frequency), and these 
regions have been mapped for both lower and upper cutoff breathers. Each region substantially 
disagreed with the theoretical prediction and, moreover, the natures of the two numerical regions 
were strongly dissimilar. The lower cutoff breathers are found to have a low-amplitude 
instability that is not predicted by the theory. The drive parameter region of the upper cutoff 
breathers contains a relatively large "finger" in which quasiperiodicity occurs, and a relatively 
large "island" in which the motion decays to rest. The corresponding instabilities are not yet 
understood, but it appears that the quasiperiodicity is a consequence of soliton shedding even 
though the shedding is not apparent except in a small region of the drive plane. This is based 
on the observation of strong quasiperiodic amplitude modulation during shedding and relatively 
weak modulation when shedding is not visible. Clearly visible soliton shedding is not expected 
to occur for a global parametric drive because a propagating breather has a spatially varying 
phase, whereas the drive is spatially constant. 
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For an end driven lattice, soliton shedding is dramatically visible. There is a drive 
amplitude threshold for the shedding to occur, and the value decreases as the linear frequency is 
approached. The phenomenon, although not yet understood, is expected to occur in any 
nonlinear wave system that possesses breathers. A possible application is the regeneration of 
fiber optic solitons which attenuate with distance. 
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APPENDIX 



A. PROGRAM COMMANDS 

1. a - Manual frequency increment, adjustable with A. 

2. A - Manual frequency increment adjustment, can be positive or negative. 

3. d - Coarse damping adjustment (decreasing), nominally .01, can be changed in 

variable declaration section of source code (the variable is Beta_increment). 

4. D - Coarse frequency adjustment (decreasing), nominally .001, can be changed in 
variable declaration section of source code (Omega increment). 

5. Ctrl D - Coarse drive amplitude adjustment (decreasing), nominally .001, can be 
changed in variable declaration section of source code (Eta_increment). 

6. e - Fine damping adjustment (decreasing), 1/10 of coarse damping increment, can 
be changed in the mainQ section of the source code under the appropriate key hit 
section. 

7. E - Fine frequency adjustment (decreasing), 1/10 of coarse frequency increment, can 
be changed in the mainQ part of the source code under the appropriate key hit 
section. 

8. Ctrl E - Fine drive amplitude adjustment (decreasing), 1/10 of coarse drive 
amplitude increment, can be changed in the mainQ section of the source code under 
the appropriate key hit section. 

9. Ctrl F - Time series record of one lattice site, prompts for a file name to store the 
series. This function works only in Text Mode (Ctrl T). The number and 
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periodicity of the samples are controlled in the eqn_motionO function section of the 
source code. Set for 80 samples, taken at the rate of one per time increment. 

10. G - Coupling adjustment (gamma), displays current value asks for new value. 

1 1 . Ctrl G - Graphics mode, displays real time lattice motion. Also refreshes the screen 
without interfering with the lattice motion. Automatically invoked when the 
program is started or restarted. 

12. Ctrl H - Total lattice energy monitor, can only be used in Text Mode. 

13. i - Zoom in. The amplification of the screen is 2” for n button pushes. The 
amplification is vertical only. 

14. Ctrl I - Increases the strobing frequency of the phaseplot option, effectively 
increasing the sampling rate. 

15. Alt K - Coupling modulation. Options are random, gradient or sine wave. 

16. Ctrl K - Amplitude kick, specified lattice site by the desired amount. Amplitude and 
lattice site input from the keyboard. 

17. Ctrl L - Pins specified lattice site. Acts like a toggle switch, reselection unpins the 
site. 

18. n - Random perturbation of the lattice. Perturbs all sites with a random coordinate 
change of < ±3% of the maximum amplitude in the lattice. 

19. o - Zoom out. The reduction of the screen is 2° for n button pushes. The reduction 
is vertical only. 

20. Ctrl O - Decreases the strobe frequency of the phaseplot option, effectively 
decreasing the sampling rate. 
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21. p - Pause program. Freezes execution of the program. Alternate method of 
obtaining a time series file. 

22. P - Phaseplot mode. Selects up to five lattice site for phase space monitoring. 

23. Ctrl Q - Exits the program. 

24. r - User initiated state save. The lattice must be stable to within 1 % (blue color) or 
the option will freeze the program until the 1 % criterion is met. The option saves 
the current lattice modulation at the upper turning point of the lattice site being 
monitored. 

25. Ctrl R - Restarts the program without going back out to DOS. Goes to the same 
screen that was displayed when the program was first started up. 

26. s - Program freeze, captures lattice modulation when pushed. Any key will continue 
the program but the lattice will be at rest. 

27. S - Monitors the stability of the specified lattice site. Works like a toggle switch. 
The lattice will change colors when the monitored site is within 1 % of its average 
amplitude for the last twenty upper turning points. 

28. Ctrl S - Turns off the drive. 

29. Ctrl T - Text Mode. Prints system parameters on screen. Total lattice energy can 
be monitored in this mode. 

30. u - Coarse damping adjustment (increasing), nominally .01, can be changed in 
variable declaration section of source code (variable is Beta_increment). 

31. Ctrl U - Coarse drive amplitude adjustment (increasing), nominally .001, can be 
changed in variable declaration section of source code (Eta_increment). 
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32. U - Coarse Frequency adjustment (increasing), nominally .001, can be changed in 
variable declaration section of source code (Omega_increment). 

33. V - Fine damping adjustment (increasing), 1/10 of coarse damping increment, can 
be changed in the mainQ section of the source code under the appropriate key hit 
section. 

34. Ctrl V - Fine drive amplitude adjustment (increasing), 1/10 of coarse frequency 
increment, can be changed in the mainQ section of the source code under the 
appropriate key hit section. 

35. V - Fine drive amplitude adjustment (increasing), 1/10 of coarse increment, can be 
changed in the mainQ part of the s urce code under the appropriate key hit section. 

36. w - Displays peak lattice amplitude, works in Graphics Mode only. 

37. Ctrl W - Turns on waterfall type display. Trends are easily noticed on this type of 
display. 

38. Ctrl X - Selects lattice site to monitor for spectrum analysis. 

39. y - Dump spectrum generated by Ctrl X. 

40. z - Sets damping to zero. 

41 . Z - Sets frequency to zero. 

42. Ctrl Z - Zeros the peak amplitude variable (variable automatically zeros when a 
parameter is changed). 
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B. PROGRAM CODE 



The QuickC program codes for the aforementioned numerical implementations (Sect III) 
are listed below. The full program listing for Lower Cutoff (CCLATL.C) is included and the 
equation of motion for Upper Cutoff (CCLATU.C), Lower Cutoff end-driven (ENDL.C) and 



Upper Cutoff end-driven (ENDU.C). 

1. Lower Cutoff, Global Drive 

I* PROGRAM LATTICE (VGA) 

VERSION 2.0 (QUICKC) 

*/ 

#define LAST UPDATE 22 JUL 1991 BY CLEON WALDEN 
I* 

THIS PROGRAM SIMULATES A GENERAL LATTICE WITH EQUATIONS OF MOTION 
THAT CAN BE SUBSTITUTED IN WHERE INDICATED. VARIOUS INTERACTIVE 
FEATURES ARE PROVIDED, WHICH ARE EXPLAINED IN THE PROGRAMMER’S 
MANUAL FOR LATTICE. ENERGY MONITORING IS ADDED TO TEXT SCREEN. 

A WATERFALL DISPLAY IS ADDED TO MONITOR SOLITON MOTION DUE TO 
MEDIUM EFFECTS. 

*1 



^include "BIOS.H" 

^include "graph. h" 

^include "stdlib.h" 

^include "CONIO.H" 

^include "MATH.H" 

^include "STDIO.H" 

^include "time.h” 

^include "direct.h" 

^include "process. h" 

^include "dos.h" 

^define getrandom(min,max) ((randQ % (int)((max)-(min)))-l-(min)-(- 1) 
#define SQR(a) ((a)*(a)) 

^define CUB(a) ((a)*(a)*(a)) 

#defme SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr 
#defme DOFOR(i,to) for(i=0;i<to;i-l- -I-) 

#defme PI 3.14159265359 
^define SCREEN_CORRECTION_F ACTOR 1.05 
^define eta_increment 0.01 
^define beta increment 0.01 
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#defme omega_increment 0.001 
#define STABILITYJNCREMENT 0.01 
#define DISPLAYJNCREMENT DINC 
#define PHASEJNCREMENT PINC 
#define MIDDLE_ELEMENT (no_pendulums/2) 

#define text_flag flags[0] 

#define graphics_flag flags[l] 

#define file_dump_flag flags(2] 

^define stability_flag flags[3] 

#define peak_flag flags[4] 

#define stable_flag flags[5] 

#define phase_flag flags[6] 

#define pause_flag flags[7] 

#define pause_flag2 flags[8] 

^define stop_flag flags[9] 

#define spectnim_flag flags[10] 

^define dump_spectrum_flag flags[ll] 

#define energy_flag flags[12] 

#define waterfall_flag flags[13] 

#define wait_flag flags[14] 

/* The dynamical variables are entered here. */ 

double coordinate[ 150],momentum[150],old_coordinate[ 150],old_momentum[ 150], 
oldold_coordinate[ 1 50] ,oldold_momentum[ 1 50] ; 

double acceleration, gamma( 150], mean_gamma,beta,omega0[150],mean_omega0, 
eta,omega,alpha,model_time,time_int, 

max_amp , mode_amp , time_ser ies_d isp [8000] ,phase_el ements [5] , 
peak_record[20],pinned_elements[150],DlNC,PINC,period, 
spectrum[4098] ,temp2 [4096] , energy , grad lent, wat_inc, 
peak_amp,amp_kick,freq_inc,delta; 

int no_pendulums,flags[ 15], counter l,phase_counter,phase_int,first_element, 
chosen_element,stability_element,disp_color,colors[5],counter2,step_size, 
spectrum_element,spectrum_counter,sample_counter,energy_counter,g,gst[8000], 
waterfall_counter,color_counter,number_clicks,p_flag,t_flag,top_flag; 

mainO { 

FILE 

int c,i,j,k,l,m,n,xx,a; 

double modulation; 

char ans[5],buffI50],buff2[50]; 

/* void display_text0,stabiiity_restart0,user_init0,display_graphics0, 
process_pause0,monitor_stability0,start_fiie_dump0, 
stop_file_dumpO,pin_elementsO,kick_elementsO,stability_checkO, 
init_phasplot0,phasplot0,compute_spectrum0,plot_spectrum0; */ 
_clearscreen(_GCLEARSCREEN); 
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user_initO; 

_setvideomode(_MRES4COLOR); 

_setcolor(l); 

sprintf(buff,''Eta %.31f Omega %.31f Gamma %.2If, eta, omega, mean_gamma); 
_registerfonts("TMSRB.FON"); 

_moveto(10,10); 

_setfont("t’tms rmn’bnS"); 

_outgtext(buff); 

XX =0; 

wat_inc=l; 
disp_color=l; 
stop_flag=0; 
counter2=0; 
energy_flag=0; 
energy _counter = 0; 
color_counter=0; 
while(stop_flag= =0) { 
energy =0; 

if(pause_flag2= = 1) { 



pause_flag2=0; 

pause_flag=0; 

process_pauseO; 

} /* if(pause...) */ 
if(kbhitO!=0) { 


/* Check for keystroke from user */ 


c=getchO; 
if(c= = 17){ 


!* 


: Quit program */ 


stop flag=l; 

} 

else if(c= = 112) { 


/* p 


: Pause program/dump state */ 


pause flag=l; 

) 

else if(c==20) { 


/* 


; Text Mode */ 



_clearscreen(_GCLEARSCREEN) ; 
text_flag = 1; 
graphics_flag=0; 
waterfall_flag= 0; 
wait_flag=0; 
spectrum_flag = 0; 

display_textO; 

} 

elseif(c==7) { I* ''G : Graphics Mode */ 

screen_graphQ; 

phase_flag=0; 

wait_flag=0; 

spectrum_flag=0; 

waterfall_flag=0; 

text_flag=0; 
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graphics_flag = 1; 

sprintf(buff,"Eta %.31f Omega %.3lf Gamma %.2ir,eta,omega,mean_gamma); 
_registerfonts('TMSRB.FON"); 

_moveto(10,10); 

_setfont("t’tms rmn’bnS"); 

_outgtext(buff); 

} 

elseif(c= = 8) { /* "'H : Monitor energy in text mode */ 

if(energy_flag= =0) energy_flag = 1; 
else energy_flag=0; 

} 

else if(c= = 16) { /* : Kick in parameter variations */ 

screen_textO; 

srand((unsigned)time(NULL)); 

printf("\nEnter 1 if you wish to vary coupling: "); 

scanf("%d",&j); 

ifO = = l){ 

printf("\nEnter 1 (random), 2 (gradient), or 3(sine) desired: "); 

scanf("%d",&k); 

if(k==2){ 

printf("\nEnter gradient (percentage over entire lattice): "); 
scanf( " % ir ,&gradient); 

DOFOR(i,no_pendulums) { 

gamma[i] = gamma[i] + mean_gamma*i*gradient/(100*no_pendulums); 

} /* DOFOR */ 

} /* if(j==2) */ 
else if(k= =3) { 

printf("\nEnter modulation amplitude (percentage): "); 
scanf(" %lf',&gradient); 

printf("\nEnter integer number of wavelengths to use: "); 
scanf("%d’’,&l); 

DOFOR(i,no_pendulums) { 

gamma[i] = gamma[i] + mean_gamma*gradient/100*(-cos(2*l*PI*i/ 
no_pendulums)); 

} 

} 

else DOFOR(i,no_pendulums) { 
j = rand0; 

gamma[i] = mean_gamma*(.95 + . 1 *j/RAND_MAX); 

} 

} 

else { 

printf("\n Varying omegao...\n"); 

printf("\nEnter 1 (random), 2 (gradient), or 3(sine) desired: "); 

scanf("%d",&k); 

if(k==2){ 

printf("\nEnter gradient (percentage over entire lattice): "); 
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scanf(" %lf',&gradient); 

DOFOR(i,no_pendulums) { 

omega0[i]=omega0[i] + mean_omega0*i*gradient/(100*no_pendulums); 

} 

} 

else if(k==3) { 

printf("\nEnter modulation amplitude (percentage): "); 
scanf(" %lf',&gradient); 

printf("\nEnter integer number of wavelengths to use: "); 
scanf("%d ",&!); 

DOFOR(i,no_pendulums) { 

omega0[i]=omega0[i] + mean_omega0*gradient/100*(-cos(2*l*PI*i/no_pendulums)); 

} 

} 

else { 

DOFOR(i,no_pendulums) { 
j = randO; 

printfC'Random number is: %d RAND_MAX is %d\n"J,RAND_MAX); 
omegaO[i] = mean_omega0*(.95 + . 1 *j/RAND_MAX); 

} 

} 

} 

screen_graphO; 

} 

else if(c= =23) { /* : Waterfall Display Mode */ 

if(waterfall_flag= =0) { 

phase_flag=0; 

graphics_flag=0; 

text_flag=0; 

spectrum_flag=0; 

waterfall_flag=l; 

initwaterfallQ; 

} 

else { 

if(wait_flag= = l) { 

wait_flag=0; 

init_waterfallO; 

} 

else waterfall_flag=0; 

} 

} 

else if(c= =80) { /* P : Phase Plot Mode */ 

if(phase_flag= =0) init_phasplotO; 
else { 

phase_flag=0; 
spectrum_flag = 0; 
graphics_flag = 1; 
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waterfall_flag = 0; 
wait_flag=0; 

screen_graphQ; 

} /* else */ 

} 

elseif(c==6) { /* : Time series to file */ 

if(file_dump_flag= =0) start_file_dumpO; 
else { 

stop_file_dumpO ; 
file_dump_flag = 0; 

} 

} 

else if(c= =47) { /* G ; Change gamma */ 

screen_textO; 

printf("\nEnter new gamma, gamma is %lf’,mean_gamma); 
scanf(" % if ' ,&mean_gamma); 

DOFOR(i,no_pendulums) { 
gamma[i] = mean_gamma; 

) 

screen_graphO; 

stability_restartO; 

} 

else if(c= = 119) { /* w ; Peak amplitude */ 

_moveto(210,10); 

_setfont("t’tms rmn’bn5"); 
sprintf(buff2,"Amp % .4if',peak_amp); 

_outgtext(buff2); 

} 

else if(c= =26) { /* : Reset peak_amp variable */ 

peak_amp=0.0; 

} 



else if(c= =24) { /* ''X ; Calculate spectrum */ 

spectrum_flag= 1; 

screen_textO; 

printf("\nEnter number of element to be analyzed: "); 
scanf(" %d'',«S:spectrum_element); 
screen_graphQ; 

} 

else if(c==21) { /* : Increase drive (coarse) */ 

save_stateO; 

eta = eta + eta_increment; 

stability _restart0; 

} 

elseif(c==4) { /* : Decrease drive (coarse)*/ 

save_stateO; 
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eta = eta-eta_increment; 
stability_restartO; 

} 

else if(c= =22) { /* : Increase drive (fine) */ 

save_stateO; 

eta = eta + eta_increment/ 1 0; 
stability_restartO; 

} 

elseif(c==5) { /* : Decrease drive (fine)*/ 

save_stateO; 

eta = eta-eta_increment/ 10; 

stability_restartO; 

} 

else if(c= =26) { /* : Turn off drive */ 

save_stateO; 

eta=0; 



stability restartO; 

} 

else if(c= =85) { /* U : Increase frequency (coarse)*/ 

save_stateO; 

time_waitO; 

time_int=time_int*200/period; 
omega = omega + omega_increment; 
period = 2 *PI/omega; 
time_int=time_int*period/200; 

stability_restart0; 

} 

else if(c= =68) { /* D ; Decrease frequency (coarse)*/ 

save_stateO; 

time_waitO; 

time_int = time_int*200/period ; 
omega = omega-omega_increment; 
period = 2 *PI/omega; 
time_int = time_int*per iod/200; 
stability_restartO; 

} 

else if(c= =86) { /* V : Increase frequency (fine)*/ 

save_stateO; 

time_waitO; 

time_int = time_int*200/period ; 
omega=omega+omega_increment/10; 
period =2*PI/omega; 
time_int = time_int*per iod/200; 
stabil ity_restartO ; 

} - 

else if(c==69) { /* E : Decrease frequency (fine)*/ 

save_stateO; 
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timewaitO; 

time_int=time_int*200/period; 
omega = omega-omegaincrement/ 1 0; 
period =2*PI/omega; 
time_int=time_int*period/200; 

stability_restartO; 

} 

else if(c==90) { /* Z : Set frequency to zero */ 

save_stateQ; 



omega =0; 

stability_restartO; 

} 

else if(c= =65) { /* : Set manual frequency increment*/ 

screen_textO; 

printfC Manual frequency increment is %lAn",freq_inc); 
printf("New increment:"); 
scanf(" % If ' ,&ffeq_inc) ; 
screen_graphO; 



} 

else if(c= =97) { /* a 

save_stateQ; 
time_waitO; 

time_int=time_int*200/period; 
omega = omega + freq_inc; 
per iod = 2 *PI/om ega; 
time_int=time_int*period/200; 

stability_restartO; 

} 

else if(c= = 1 17) { /* u 

save_stateO; 
beta = beta + beta_increment; 
stability_restartO; 

} 

else if(c= = 100) { /* d 

save_stateQ; 

beta = beta-beta_increment; 

stability_restartO; 

} 

else if(c= = 118) { /* v 

save_stateO; 

beta = beta + betaincrement/ 10; 
stability restartO; 

} 

else if(c= = 121) { /* y 

if(spectrum_flag= = 1) dump_spectrum_flag = 1; 

} 



Manual Frequency increment*/ 



Increase damping (coarse)*/ 



Decrease damping (coarse)*/ 



Increase damping (fine)*/ 



Dump spectrum */ 



else if(c= = 101) { /* e 



Decrease damping (fine)*/ 
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savestateQ; 

beta = beta-beta_increment/ 1 0; 
stability_restartO; 

} 

else if(c= = 122) { /* z : Turn off damping */ 

save_stateO; 



beta=0; 

stability_restartO; 

} 

else if(c= = 12) { /* 

pin_elementsO; 

stability_restartO; 

} 

else if(c= = ll) { /* "K 

kick_elementsO; 

stability_restartO; 

} 

else if(c= = 105) { /* i 

DINC=DINC/2; 
P1NC=P1NC*2; 
watjnc = wat_inc*2; 

} 

else if(c= = 1 1 1) { /* o 

DINC = DINC*2; 
PlNC=PlNC/2; 
wat_inc = wat_inc/2; 



: Pin elements */ 



: Kick elements */ 



: Zoom in */ 



: Zoom out */ 



} 

else if(c= = 1 10) { /* n : Perturb system */ 

I* Maximum +/-3% of peak value */ 
linear_kickO; 
screen_graphO; 
stability_restartO; 



} 

elseif(c= = 15) { /* "O 

phase_int = phase_int- 1 ; 

} 

elseif(c==9) { /* "I 

phase_int=phase_int+ 1 ; 

} 

elseif(c= = 18) { /* 

screen_textO; 
userJnitO; 

sprintf(buff,"Eta %.31f Omega %.31f Gamma %.21f', eta, omega, mean_gamma); 
_registerfonts(’TMSRB.FON"); 

_moveto(10,10); 

_setfont("t’tms rmn’bn5"); 

_outgtext(buff); 



: Decrease strobe freq */ 
Increase strobe freq */ 

; Restart program *t 
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if(phase_flag= =0) _setvideomode(_MRES4COLOR); 
else _setvideomode(_VRES16COLOR); 

} 

else if(c= = 115) { /* s : Stop everything as is */ 

DOFOR(i,no_pendulums) { 

coordinate[i]=0; 

momentum[i]=0; 

} I* DOFOR */ 
stability_restartO; 

} 

else if(c= = 83) { /* S : Monitor stability */ 

mo nitorstabil ity 0; 

sprintf(buff,"Eta %.31f Omega %.31f Gamma %.21f',eta,omega,mean_gamma); 
_registerfonts("TMSRB.FON"); 

_moveto(10,10); 

_setfont("t’tms rmn’bnS"); 

_outgtext(bufO; 

} 

else if(c= = 1 14) { /* r : User initiated state save *! 

save_stateuO; 

} 

} /* if(kbhit...) */ 

if(wait_flag= =0) { 
eqn_motion0; 

} I* if(wait_flag */ 

} I* while(stop_flag... *! 

_setvideomode(_DEFAULTMODE); 
printffPROGRAM COMPLETE AT %lf',model_time); 

} I* MAINO */ 

user_initO { 

FILE *fq; 
int c,i,j,k,l,m; 

char answer[l],ans[l],answ[l],filename[30]; 

printfC GENERALIZED LATTICE MODEL PROGRAM W/ VARIABLE PARAMETERSVn"); 

printffVersion 2.1X486 (MS) \n"); 

printfC'Last updated 22 JUL 1991\n"); 

printfC'Variant notes: Default IC is AM\n"); 

printf("OmegaO is set equal to one for all cases!\n"); 

printf(" Rotating phase plane is used.\n"); 

printf("Real time FFT function is added... \n"); 

printf("Energy monitoring available via in text modeVn"); 

printfC'Set top_flag to get output files for theory comparison.\n"); 

printf(”Manual Frequency increment is set to 1/20 omega inc.\n”); 

/* Sets flag to output additional information for theory comparison. 

Adds the information to the end of the normal output file*/ 



64 



printfC'Set top_flag? "); 
scanf("%s",ans); 

if((ans[0]==89)l |(ans[0] = = 121)) { 

top_flag=l;} 

else top_flag=0; 

printf("\nDo you want to use a file for initial conditions (Y/N)? "); 
scanf("%s", answer); 
sampl e_cou nter = 0; 

if((answer[0] = = 89) j 1 (answer[0] = = 121)) { 

printf("\n01d file?"); 

scanf("%s",answ); 

printf("\nEnter name of file to be read: "); 
scanf(" %s", filename); 

if((fq=fopen(filename,"r"))! =NULL) { 
mean_gamma=0; 

fscanf(fq,"%d\n",&no_pendulums); 
fscanf(fq,"%lAn",&alpha); 
fscanf(fq," %lAn",&beta); 
fscanf(fq," %lAn" ,&eta); 
fscanf(fq," %lAn",&omega); 

DOFOR(i,no_pendulums) { 
fscanf(fq,"%lf %lf %lf %lAn", 
&omegaO[i],&gamma[i],&coordinate[i],&momentum[i]); 
mean_gamma = meangamma + gamma[ i] ; 

} 

mean_gamma = mean_gamma/no_pendulums ; 
fclose(fq); 

} 

else printfC'Can’t open file requested."); 

} /* if((ans... */ 
else { 

printf("\nEnter number of pendulums to use: "); 
scanf(" %d",&no_pendulums); 
printf("\nEnter mode amplitude: "); 
scanf(" % If ' ,&mode_amp); 
printf("\nEnter modulation amplitude: "); 
scanf(" % If ,&max_amp); 

printf("\nEnter nonlinear coefficient alpha (+/- 1 ONLY): "); 
scanf(" % If ' ,&alpha); 

DOFOR(k,no_pendulums) { 

if(alpha< 0) coordinate[k] =pow((-l),k)*(mode_amp 

+max_amp*sin(2*PI*k/no_pendulums)); 

else coordinate[k] =mode_amp+max_amp*sin(2*PI*k/no_pendulums); 

momentum[k]=0; 

pinned_elements[k] =0; 

} /* DOFOR */ 

printf("\nEnter coupling coefficient gamma: "); 
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scanf(" %lf’,«&meaii_ganuna); 

DOFOR(i,no_pendulums) { 
gamma[i] = mean_gamma; 
omegaO(i] = 1; 

} 

printf("\nEnter drive amplitude eta: "); 
scanf(" %lf',«&eta); 

printf("\nEnter drive frequency omega: "); 
scanf(" %lf',&omega); 

printf("\nEnter dissipation constant beta: "); 
scanf( " % If ' ,&beta) ; 

} /* else */ 

printf("\nEnter number of steps per response cycle: "); 
scanf(" %d",&step_size); 
if ((step_size/2) != 0) g=.5; 
else g=0; 

time_int=2*PI/(step_size*omega); 

DOFOR(i,15) flags[i]=0; 

DOFOR(i,4096) spectrum[i] = 0; 

DOFOR(i,150) { 
oldold_coordinate[i] = 0.0; 
old_coordinate[i] =0.0; 

} 

screen_graph0; 
number_clicks=0; 
mean_omega0= 1 .0; 
peak_amp=0; 
spectrum_counter = 0; 

/* Phase correction to start files in the proper phase*/ 

if((answ(0] = = 89) j | (answ(0] = = 121)) delta =asin((omega*beta)/eta)/2; 

else delta = (asin((omega*beta)/eta)/2)-PI; 

model_time=0.0; 

freq_inc= omega_increment/20; 

DINC = .02; /* CHANGE THIS TO CHANGE SCALE OF DISPLAY */ 

PINC=2; /* CHANGE THIS TO CHANGE SCALE OF PHASE PLOT */ 
} /* USERJNIT */ 



display_graphics0 { 

/* Displays the lattice sites. The screen updates sequentially through the 
points after each round of motion calculations are complete. */ 
int c,ij,k,l,m,n; 

if(no_pendulums < =40) { 
first_element=0; 

DOFOR(k,no_pendulums) { 
n=0; 

if((stability_element= =k)&&(stability_flag= = 1)) n= 1; 
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if((pinned_elements[k] = = l)&&(disp_color= = 1)) n=2; 
if((pinned_elements[k] = = l)&&(disp_color= =2)) n=-l; 
l=old_coordinate[k]/DISPLAY_INCREMENT; 

_setcolor(0); 

_setpixel((160-5*MIDDLE_ELEMENT + 5*k), (100+1)); 
l=coordinate[k]/DISPLAY_INCREMENT; 

_setcolor(d isp_coior + n); 

_setpixel((160-5*MlDDLE_ELEMENT + 5 *k), (100+1)); 

} /* DOFOR */ 

} /* IF */ 
else { 

1 = (no_pendulums-40)/2; 
first_element=l; 

DOFOR(k,40) { 

m = old_coordinate[I + k]/DISPLA YJNCREMENT; 
n=0; 

if((stability_element= =(l + k))&&(stability_flag= = 1)) n= 1; 
if((pinned_elements[k]= = !)&&(disp_color= = 1)) n=2; 
if((pinned_elements[k] = = !)&&(disp_color= =2)) n=-l; 
_setcolor(0); 

_setpixel((60 + 5*k),(100+m)); 
m=coordinate[l + k]/DISPLA YJNCREMENT; 
_setcolor(disp_color + n); 

_setpixei((60 + 5*k),(l(X)+m)); 

} /* DOFOR */ 

} /* ELSE */ 

} I* DISPLAY_GRAPH1CS *! 



display JextQ { 

/* The coordinate and other information below is a snapshot of what the 
lattice is like at the particular instant in time when the key was pushed.*/ 
char message[80]; 
int c,j,k,l,m; 

_setvideomode(_DEFAULTMODE); 

printf("Time is : % If, model jime); 

printf(" System parameters are: \n"); 

printf(" Gamma %lf',mean_gamma); 

printf(" Eta %lf',eta); 

printfC Omega %lf\n",omega); 

printf(" Beta %lf',beta); 

printf(" Alpha %lAn", alpha); 

printf("\nThere are %d elements in the system", nojpendulums); 

printf("\n\nPress any key to continue..."); 

c=getcharO; 

while(kbhitO= = 0 ) ; 

printfC'Element Position Velocity Element Position Velocity\n"); 
DOFOR(j,20) { 
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%lAn", 



printfC %d %lf %lf %d %lf 

j,coordinate[first_element+j],momentum[first_element+j], 

(j + 20),coordinate[first_element+ j + 20], 
momentum[first_element + j + 20]); 

} /* DOFOR */ 

printf("\n\nPress for graphics, to quit."); 
c=getcharO; 
while(kbhitO= =0); 
while(kbhitO= =0); 

} /* display_text */ 

process_pauseO { 

/* Old method of saving files. Not sure if it has other uses.*/ 
int c,i,j,k,l,mode; 

FILE *fr; 

char filename[30], message[80]; 

_ciearscreen(_GCLEARSCREEN) ; 

_setvideomode(_DEF AU LTMODE) ; 
printfC'Do you wish to save this state? "); 
if((c=getcheO)= = 121) { 
printfC \iiEnter name of file to be written: “); 
scanfC %s", filename); 

if((fr= fopen(filename,"w"))! = NULL) { 
fprintf(fr, " %d\n" ,no_pendulums); 
fprintf(lf , " %1 An" ,alpha); 
f])rintf(fr," %IAn" ,beta); 
fprintf(fr,"%lAn",eta); 
f])rintf(ff , " % 1 An " , omega); 

DOFOR(i,no_pendulums) 
fprintf(fr,"%lf %lf %lf %lAn", 
omegaO[i],gamma[i],old_coordinate[i],momentum[i]); 

fclose(fr); 

} /* ifO */ 

else printf("Failed to open %s\n", filename); 

} /* if (0) */ 

time_int = time_int*200/per iod ; 
printfC \nEnter new time multiple (old multiple is %lf): ",time_int); 
scanfC %lf',&time_int); 
timeint = time_int*per iod/2(K); 
if(phase_flag= =0) _setvideomode(_MRES4COLOR); 
else _setvideomode(_VRES16COLOR); 

} /* process_pause */ 

start_file_dump0 { 

/* Sets up a file dump of the coordinates of the selected element. Set up 
to record every time through for 30 times */ 
_clearscreen(_GCLEARSCREEN); 
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_setvideomode(_DEFAULTMODE); 
printf("Enter number of element to be monitored: "); 
scanf(" % d " ,&chosen_element); 
if(chosen_element > (no_pendulums-l)) 
printf("Out of range. No file dump"); 
else file_dump_flag = 1; 

_clearscreen(_GCLEARSCREEN); 
if(phase_flag= =0) 

_setvideomode(_MRES4COLOR); 

else 

_setvideomode(_VRES 16COLOR); 
counter 1 =0; 

} /* start_file_dump */ 

stop_file_dumpO { 

/* Dumps the collected data to a file. */ 
char filename[30]; 
int c; 

FILE *fs; 

_cl earscreen(_GCLE ARS CREEN) ; 
_setvideomode(_DEFAULTMODE); 
printf("\nEnter name of file to be written: "); 
scanf("%s", filename); 
if((fs=fopen(filename,"w"))! = NULL) { 

DOFOR(c,8000) fprintf(fs,"%lf %dfin",time_series_disp[c],gst(c]); 

_clearscreen(_GCLEARSCREEN); 

counter 1 =0; 

if(phase_flag= =0) 

_setvid eomod e(_MRES4COLOR) ; 
else _setvideomode(_VRES16COLOR); 

} /* if(0) */ 

} /* stop_file_dump */ 

monitor stabilityO { 

/* Selects the lattice site to be monitored */ 
if(stability_flag= =0) { 
stability_flag=l; 
disp_color=2; 

_setvideomode(_DEFAULTMODE); 
printf("Enter number of element to be monitored: "); 
scanf(" %d",&stability_element); 
stability_restart Oi 
} 

else { 

stability_flag=0; 

disp_color=l; 

} 
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if(phase_flag= =0) 

_setvideomode(_MRES4COLOR); 

else 

_setvideomode(_VRES 16COLOR); 

} 

stability_checkO { 

/* The criterion is set for 1 % . Appears to check to see if 20 consecutive 
peaks are within the criterion of each other. */ 
int i,k; 
k=0; 

DOFOR(i,19) { 
stable_flag= 1; 

if((peak_record[i+ 1] > (peak_record[i]*(H-STABILITY_INCREMENT))) | 
(peak_record[i+l]<(peak_record[i]*(l - STABILITYJNCREMENT )))) { 
stable_flag=0; 
break; 

} /* if */ 

} I* DOFOR */ 
if(stable_flag= = 1) { 
disp_color= 1; 

} 

} 

pin_elementsO { 

/* Pins an element at rest */ 
int ans,i,check; 

_setvideomode(_DEFAULTMODE); 
check =0; 

printf("\nEnter number of element to be pinned: "); 
scanf("%d",&ans); 

if((ans < 0) ] (ans > (no_pendulums- 1 ))) { 
check= 1; 

} 

if(check==0) { 
if(pinned_elements[ans]= =0) { 
pinned_elements[ans] = 1 ; 
coordinate[ans] = 0; 
momentum[ans] =0; 

} 

else pinned_elements[ans]=0; 

} /* if *! 

if(phase_flag= =0) 

_setvideomode(_MRES4COLOR) ; 
else 

_setvideomode(_VRES 16COLOR); 

} 
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kick_elementsO { 

I* Equivalent to banging an element with a hammer. No control over when 
the hit occurs in the cycle. */ 
int ans, check; 
double amount; 

_setvideomode(_DEF AULTMODE); 
check =0; 

printf("\nEnter number of element to kick: "); 
scanf("%d",&ans); 

if((ans<0)|(ans>(no_pendulums-l))) { 
check = 1; 

} 

if(check==0) { 

printf("\n Enter amount to kick: "); 

scanf(" % If ' , (Seamount) ; 

coord inate[ans] = coordinate[ans] + amount; 

} /* if */ 

if(phase_flag= =0) 

_setvideomode(_MRES4COLOR); 

else 

_setvideomode(_VRES 16COLOR); 

} 

stability_restartO { 

/* Changes the color to show restart and zeros the stability check matrix */ 
int i; 

char buffi [50]; 
disp_color=2; 
stable_flag=0; 

DOFOR(i,20) peak_record[i]=0; 
peakrecordf 15] = 10; 
peak_amp=0; 

sprintf(buffl,"Eta %.31f Omega %.31f Gamma %.21f',eta,omega,mean_gamma); 
_moveto(10,10); 

_setfont("t’tms rnm’bn5"); 

_outgtext(buffl); 

r 

save_stateO { 

/* Automatic state saver. It’s supposed to keep the program from crashing 
when a lattice flys off numerically. */ 
int i; 

char filename[] = ’’savedsta"; 

FILE *fl; 

I* change the number below to change strokes between autosaves */ 
if(number_clicks = = 100) { 
if(chdir("E:\\MODEL\\CLEON") != 0) { 
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screen_textO; 

printfC'failed to change to E:\\MODEL\\CLEON.\n"); 
stability restartO; } 

p _flag=0; 
while(p_flag= =0) { 
eqn_motionO; 

} 

if((fl=fopen(filename,"w")) != NULL) { 

fprintf(fl," %d\n",no J5endulums); 
fprintf(fl,"%mn",alpha); 
fprintf(fl,"%lAn",beta); 
fprintf(fl,"%mn",eta); 
f^rintf(fl ,"% 1 An" , omega); 

DOFOR(i,no_pendulums) 
fprintf(fl,"%lf %If %lf %lAn", 
omegaO[i],gamma[i],old_coordinate[i],old_momentumIi]); 
fclose(fl); 

number_clicks=0; } 
else { screen_textO; 

printfC'Failed to open %s\n", filename); 
stability_restart; } } 
if(chdir("E:\\MODEL") != 0) { 
screen_textO; 

printfC'failed to change back to E:\\MODEL.\n"); 

stability_restartO; } 

else + +number_clicks; 
screen_graphO; } 

/* save_state */ 
save_stateuO { 

/* User initiated state save. Normal way to save lattice data files. Saves 
data for the lattice at a turning point (within one time step). */ 
int i; 

char filename[30]; 

FILE *Ip; 
screen_textO; 

printf("\nEnter name of file to save state in: "); 
scanf(" %s",filename); 
p_flag=0; 
while(p_flag= =0) { 
eqn_motion0; 

} 

printf("\n%lf % If ,model_time, delta); 
if((Ip=fopen(filename,"w")) != NULL) { 
fprintf(fp, " % d\n" ,no_pendulums); 
fprintf(fp," %lAn",alpha); 

I^r intf(f^, " % 1 An " ,beta); 
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fprintf(fp,"%lf\n",eta): 
fprintf(fi), "% 1 An" , omega); 

DOFOR(i,no_pendulums) 

fprintf(fp,"%lf %If %lf %lAn", 
omegaO[i],gamma(i],old_coordinate[i],old_momentum[i]); 
if(top_flag= = 1){ 

fprintf(fp,"%lf %lf %lf %lAn",oldold_coordinate[stability_element], 
old_coordinate[stability_element],coordinate[stabiIity_element],time_int); 

fprintf(fp,"%lf %lf %lf',oldold_momentum[stability_element], 
old_momentum[stability_element],momentum[stability_element]);} 
fclose(fp); } 

else printfO'Failed to open %s\n", filename); 
pause_flag=0; 

printf("\nEnter number of steps per response cycle: "); 
scanf(" %d" ,&step_size); 
if ((step_size%2) != 0) g=.5; 
else g=0; 

timejnt = 2 *PI/(step_size*omega); 
screen_graphO: } 

/* save_stateu */ 

retrieve_stateO { 

/* Retrieves the auto saved file when the lattice crashes. */ 
int i,c; 

char filenameG = "savedsta"; 

FILE *fg; 
screen_textO; 

printf("\nLATTICE CRASH! !!!\n"); 

printf("Do you want to retrieve the latest stored state(Y/N).\n"); 
printf(" Eta %lf Omega %lAn",eta, omega); 
if((c=getche0)= = 121) { 
if(chdir("E:\\MODEL\\CLEON") != 0) { 
printf("Failed to change directories. \n"); 

} 

} 

else { printf("\nEnter name of file to retrieve.: "); 

scanf("%s",filename); } 
if((fg=fopen(filename,"r"))! =NULL) { 
mean_gamma=0; 

fscanf(fg," %d\n",&no_pendulums); 
fscanf(fg, " % lAn" ,«S:alpha); 
fscanf(fg,"%lAn",&beta); 
fscanf(fg,"%lAn",&eta); 
fscanf(fg," %lAn",&x)mega); 

DOFOR(i,no_pendulums) { 
fscanf(fg,"%lf %lf %lf %lAn", 
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&omegaO[i],&gamma[i],&coordinate[i],&momentum[i]); 
mean_gamma = mean_gamma + gamma[i] ; 

} 

fclose(fg); 

counter2=0; 

delta = atan((omega*beta)/sqrt(SQR(eta)-(SQR(omega)*SQR(beta))))/2; 
model_time=0.0; 

} 

else printf("Can’t open file requested."); 
time_int=time_int*200/period; 

printf("\nEnter new time multiple (old multiple is %lf): ",time_int); 

scanf(" %lf',&time_int); 

time_int=time_int*period/200; 

if(chdir("E:\\MODEL") != 0) { 

screen_textO; 

printf("failed to change back to E;\\MODEL.\n"); } 
screen_graphO; 

} 

I* retrieve_state */ 
screen_textQ { 

_clearscreen(_GCLEARSCREEN); 
text_flag= 1; 
graphics_flag=0; 
spectrum_flag=0; 

_setvideomode(_DEFAULTMODE); } 

/* screen_text *! 

screen_graphO { 

_clearscreen(_GCLEARSCREEN); 

_setvideomode(_MRES4COLOR); 

_setcolor(l); 
phase_flag=0; 
spectrum_flag=0; 
text_flag=0; 
graphics_flag=l; } 

/* screen_graph *! 

Iinear_kick0 { 

/* Lattice perturbation up to a maximum of 3%. Sort of random as to sign of 
kick direction for each site. Uses the % of the maximum amplitude site. */ 
int aj,i; 
double k,b; 

struct dostime_t time2; 

_dos_gettime(&time2); 
srand(time2 .hsecond); 

DOFOR(i,no_pendulums) { 
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b=randO; 

/* change 33.33 to change %, now 3% */ 
k=b/(33.33 * 32767); 
amp_kick= k*peak_amp; 
if(b > 16384) { 

coordinate[i) = coordinate[i] +amp_kick; 

} 

else coordinate[i] = coordinate[i] - amp_kick; 

} 

peak_amp=0.0; 

} /* linear_kick */ 

time_waitO { 

/* part of a routine to get the lattice at the turning point. */ 
t_flag=0; 
while(t_flag= =0){ 
eqn_motionQ; 

} 

}/* time_wait */ 

init_phasplotO { 

/* Sets up the phase plot to look at the phase plane of sites. */ 
int c,i,j,k,l; 

char ans[5],message[40],txt[3]; 
float dummy; 

spectrum_flag=0; 

_setvideomode(_DEFAULTMODE); 

DOFOR(i,5) phase_elements[i]=0; 
printf(”\nDo you want Poincare sections? "); 
if((c=getchO)= = 121) phase_flag=2; 
else phase_flag = 1; 

printf("\nWhich elements do you wish to monitor (999 to finish): "); 

i=-i; 

k=0; 

while((i+ + ! = 999)&&(k<5)) { 
scanf(”%d",&i); 
phase_elements[k + + ] = i + 1 ; 

} /* while */ 

DOFOR(i,5) { 

if(phase_elements[i] = = 1000) { 
phase_elements[i] =0; 

} 

} 

_setvideomode(_VRES 16COLOR); 
_clearscreen(_GCLEARSCREEN); 
graphics_flag=0; 
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_settextcolor(7); 
colors[0]=2; 
colors[l]=3; 
colors[2]=9; 
colors[3] = 14; 
colors[4] = 13; 

DOFOR(i,120) { 

_setpixel(318,4*i); 

} 

DOFOR(i,160) { 

_setpixel(4*i,238); 

} 

D0F0R(i,6) { 

_setpixel(3 19,1 0*SCREEN_CORRECTlON_FACTOR + 1 
+ (80/SCREEN_CORRECTION_FACTOR)*(i + 1)); 

} 

D0F0R(i,8) { 

_setpixel(80*(i+ 1),239); 

} 

_moveto(480,450) ; 

D0F0R(i,5) { 
if((phase_elements[i])!=0) { 
c=phase_elements[i]-l; 

/* itoa(c,txt,10); 

_outtext(txt); */ 
printf("%d ”,c); 

_setcolor(colors[i]); 

DOFOR0,4) { 

D0F0R(k,8) { 

_setpixel(550+ 10*i + k, 10+j); 

} 

} 

} 

} /* DOFOR */ 
printfC' LAT2X486.C"); 

printf(”\nGamma: %lf Eta; %lf Beta: %lf',mean_gamma,eta,beta); 

printf("\nOmega: %lf Alpha: %lf Time interval; %lf, 

omega, alpha, timejnt); 

dummy= l/(time_int*omega); 

phasejnt = dummy/ 1 ; 

phase_counter=0; 

} 

phasplotO { 

/* Displays the phase for chosen elements. */ 
int i,j,k,l; 

double speed, position,fast_coordinate,fast_momentum; 
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if(phase_flag= = 1) { 

D0F0R(i,5) { 

if((j=phase_elements[i]-l)! = -l) { 

speed = PH ASE_INCREMENT*momentum(j]; 

_setcolor(colors[i]); 

position=PHASE_INCREMENT*coordinate[j]; 
fast_coordinate = (position*cos(omega*model_time) 

-speed *sin(omega*model_time)/omega) 

/(3*SCREEN_C0RRECTI0N_FACT0R); 

fast_momentuin=(position*sin(omega*model_time) 

-l-speed*cos(omega*model_time)/omega)/4; 

_setpixel((320 + 320*fast_coordinate),(240 -I- 240*fast_momentum)); 

} /* if(0) */ 

} /* DOFOR *! 

} I* if */ 

else if(phase_flag= =2) { 
if(phase_counter= =phase_int) { 
phase_counter=0; 

DOFOR(i,5) { 

if((j=phase_elements[i]-l)!=-l) { 

speed = PH ASE_INCREMENT*momentum[j]; 

_setcolor(colors[i]); 

position = PH ASE_INCREMENT*coordinate[j] ; 
fast_coordinate = (position*cos(omega*model_time) 
-speed*sin(omega*model_time)/omega); 
fast_momentum = (position*sin(omega*model_tinie) 

+ speed *cos(omega*model_time)/omega)/4; 
_setpixel((320-l-320*fast_coordinate),(240-t-240*fast_momentum)); 

} /* if(0) */ 

} /* DOFOR */ 

} /* if */ 

} I* else ifO */ 

} /* phasplot */ 

/* FFT Routines for lattice programs */ 

compute_spectrumO{ /* Uses algorithm from p. 41 1 of Numerical Recipes in C */ 

int n,mmax,mj,istep,i,nn,temp; 

double wtemp,wr,wpr,wpi,wi,theta,tempr, tempi; 

nn=2048; 

n = nn< < 1; 

j=i; 

DOFOR(i,2048) { 
nn = (i&1024)/1024; 
nn=nn-l-(i«&512)/256; 
im = nn + (i&256)/64; 
nn=nn+(i&128)/16; 
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nn=nn+(i&64)/4; 

nn=nn+(i&32); 

nn = nn + (i& 1 6)*4 ; 

nn=nn+(i&8)*16; 

nn = nn + (i&4)*64; 

nn = nn + (i&2)*256; 

nn = nn + (i& 1 )* 1 024; 

temp2[2 *nn] = spectrum(2 *i] ; 

temp2[2*nn+ l]=spectrum[2*i + 1]; 

} 

DOFOR(i,4096) spectrum[i]=temp2[i]; 

mmax=2; 

while(n> mmax) { 

istep=2*mmax; 

theta = 2 *PI/mmax ; 

wtemp = sin(0.5 *theta); 

wpr = -2.0*wtemp*wtemp; 

wpi=sin(theta); 

wr = 1.0; 

wi=0.0; 

for(m=0;m<(mmax-l);nH- =2) { 
for(i=m;i< =n;i+=istep) { 
j = i+mmax; 

tempr = wr*spectrum[j] -wi*spectrum[j + 1 ] ; 
tempi = wr*spectrum[j + 1 ] + wi*spectrum[j] ; 
spectrum(j] = spectrum[i]-tempr; 
spectrum(j + l]=spectrum[i+ l]-tempi; 
spectrum[i]+ =tempr; 
spectrum[i+ 1]+ =tempi; 

} 

wr = (wtemp = wr)*wpr-wi*wpi + wr; 
w i = wi* wpr + wtemp*wpi + wi ; 

} 

mmax = istep; 

} 

} 

plot_spectrumO { 
int ij; 

double freq, magnitude, mag, max; 
int c; 

FILE *fp; 

if(dump_spectrum_flag= = l) { 
printf("\nDumping spectrum now"); 
fp=fopen("spectrum.out","w"); 
for(i=l;i<2048;i++) { 
if(time_int! =0) fireq = i/(2048*time_int); 
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else printf("Time_int was zero!"); 
magnitude=sqrt(spectrum[(2*i)]*spectrum[(2‘*‘i)] 

+ spectrum[(2*i+ l)]*spectrum[(2*i+ 1)]); 
fprintf(fp,"%If %lf \n",fTeq,magnitude); 

} 

fclose(fp); 

dump_spectrum_flag = 0; 

} 

_clearscreen(_GCLEARSCREEN); 

max=0; 

for(i=l;i<2049;i++) { /* this loop calculates the max spectral component*/ 
magnitude=sqrt(spectrum[(2*(i))]*spectrum[(2*(i))] 

+spectrum[(2*(i)+ l)]*spectrum[(2*(i)+ 1)]); 
if(magnitude > max) max = magnitude; 

} 

_setcolor(3); 

DOFOR(i,640) _setpixel(i,460); 

DOFOR(i,l 1) _setpixel(51*i,461),_setpixel(51*i,462); 

printfC'Spectrum for element %d LATTIC2X.C",spectrum_element); 

printf("\nAlpha: %lfBeta: % If Gamma: %lfEta: % If Omega: %lf', 

alpha,beta,mean_gamma,eta,omega); 

printf("\nMax amp = %lf ”,max); 

freq=512/(2048*time_int); 

printf("Max freq shown: %lf Time interval: %lf',freq,time_int); 
_setcolor(l); 

DOFOR(i,512) { 
j=0; 

magnitude=sqrt(spectrum[2*i]*spectrum[2*i] 

+ spectrum[2*i+ l]*spectrum[2*i+ 1]); 
mag = magnitude*400/max; 
while(j + 4- < mag) { 

_setpixel(i,460-j); 

} /* while */ 

} /* DOFOR */ 

} 

init_waterfallO { 

/* Sets up waterfall type display to see trends in latice motion. Good for 
observing kink or breather drift. */ 
int c,i,j,k,l,m; 
char ans[5]; 

_setvideomode(_VRES 1 6COLOR); 

_setcolor(l); 

DOFOR(i,600) _setpixel(20 + i,20); 

DOFOR(i,450) _setpixel(20,i4-20); 

DOFOR(i,61) { 

if((i = = 10) I (i = = 20) |(i = = 32) |(i = = 42) I (i = = 52)) _setcolor(2); 
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D0F0R(j,113) { 

_setpixel(20+ 10*i,20+4*j); 

} 

_setcolor(l); 

} 

_setcolor(2); 

DOFOR(i,ll) { 

DOFORO,160) { 

_setpixel(20+4*j,16+40*(i+ 1)); 

} 

} 

waterfall_counter = 0; 
color_counter=2; 

} 

disp_waterfallO { 
int 

if(waterfall_counter+ + < 1 10) { 
if(color_counter+ + = = 14) { 
color_counter=2; 

_setcolor(color_counter); 

} 

else { 

_setcolor(color_counter); 

} 

DOFOR(i,no_pendulums) { 
if(waterfall_counter<55) { 

_setpixel(20 + 2 *i ,24 + 8 *waterfall_counter-7*wat_inc*(coord inate[i])) ; 

} 

else { 

_setpixel (340+2*1,24+ 8*(waterfall_counter-55)-7*wat_inc*(coordinate[i])); 

} 

} 

} 

else { 

wait_flag = 1; 

} 

} 

eqn_motion0 { 

/* The following loop calculates a round of velocities for one 
iteration of the model clock. This is where the equations 
of motion need to be located, if one wishes to change the 
model’s physics! */ 
int i,j,k; 

DOF OR(i, no_pendulums) 
if(coordinate[i] > 100) { 
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retrieve_stateO; 

} 

DOFOR(i,no_pendulums) { 
if(i==0){ 

acceleration = gamma[i] *(coordinate[i + 1] + coordinate[no_pendulums- 1 ] 
-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])+2*eta*cos(2*omega*model_time-delta))*coordinate[i] 
+ alpha*CUB(coord inate[i]) ; 

} 

else if(i= = (no_pendulums-l)) { 

acceleration = gamma[i]*(coordinate[i-l] + coordinate[0] 

-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])+2*eta*cos(2*omega*model_time-delta))*coordinate[i] 

+alpha*CUB(coordinate[i]); 

} 

else { 

acceleration = gamma[i] *(coordinate[i + 1 ] + coordinate[i- 1 ] 
-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])+2*eta*cos(2*omega*model_time-delta))*coordinate[i] 
+ alpha*CUB(coordinate[i]); 

} 

oldold_momentum[i]=old_momentum[i]; 

old_momentum[i] =momentum[i]; 

momentum[i] =momentum[i] +acceleration*time_int; 

} /* DOFOR */ 

DOFOR(i,no_pendulums) { 
oldold_coord inate[ i] = old_coord inate[ i] ; 
old_coordinate[i] = coordinate[i]; 
if(pinned_elements[i] = =0) { 
coordinate[i] =coordinate[i] + momentum[i]*time_int; 

} 

/* Records data for file dump */ 

if((i= =chosen_element)&&(file_dump_flag= = 1)) { 

I* if(counter2+ + = =30) { */ 
time_series_disp[counterl + +]=coordinate[i]; 
gst[counterl + + ]=g; 

/* counter2=0; 

}*/ 

if(counterl = =8000) stop_file_dumpO; 

} /* if((i... */ 

I* Stores energy information */ 

if(energy_flag= = l) { 

energy = energy + 0. 5 *(S QR(momentu m[ i] ) 

+ gamma[i]*SQR((coordinate[i+ 1] -coordinate! i])) 
-l-SQR(coordinate[i]))-alpha/4*pow(coordinate[i],4); 
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} /* if(energy...) */ 



/* Stores FFT information */ 
if((spectrum_flag= = l)&&(spectrum_element= =i) 
&&((sample_counter+ + )= = !)) { 
if(spectrum_counter < 2048) { 
spectrum[2*spectrum_counter] = coordinate! i]; 
spectrum[2*spectrum_counter + 1] =0; 
spectrum_counter + + ; 
sample_counter = 0; 

} 

else { 

compute_spectrumO; 

_clearscreen(_GCLEARSCREEN); 

_setvideomode(_VRES 16COLOR); 

plot_spectrumO; 

phase_flag=0; 

text_flag=0; 

graph ics_flag=0; 

DOFOR(j,4096) spectrum[j]=0; 
spectrum_counter = 0; 
sampl e_counter =0; 

} /* else */ 

} /* if ((spec... */ 

/* Records the peaks for the stability checker */ 

if(((i= =stability_element)&&(stability_flag= = l)&&(stable_flag= =0))| 
((i = =(no_pendulums-l))&&(waterfall_flag= = 1))) { 
if((coordinate[i] >0)&&(coordinate[i] <old_coordinate[i]) 
&&(peak_flag==0)) { 
peak_flag=l; 

if(pause_flag= = 1) pause_flag2 = 1 ; 

DOFOR(k,19) peak_record[k]=peak_record[k+ 1]; 
peak_record[ 19] =coordinate[i]; 
if(stability_flag= = 1) stability_checkO; 
if(waterfall_flag= = 1) disp_waterfallO; 

} 

else if(coordinate[i] < 0) peak_flag=0; 

} 

/* Records peak amplitude for display and linear kick */ 
if((coordinate[i] > 0)&&(coordinate[i] < old_coordinate[i])) { 
if(old_coordinate[i] > peak_amp){ 
peak_amp = old_coordinate[i]; 

} 

} 

if((i= =stability_element)&&(stability_flag= = l)&&(stable_flag= = 1)){ 



82 



if(((coordinate[i]-old_coor(iinate[i]) > = 0) && 

((old_coordinate[i]-oldold_coordinate[i]) < = 0) && (coordinate[i] < 0)){ 
p_flag=l; 

} 

} 

} /* DOFOR */ 

/* Increments time and resets it at 2Pi. Causes problems in eqns of 

motion if not reset */ 

model_time = mod el_time + 1 ime_int ; 

if((g+l) = = (step_size/2)){ 

model_time=0.0; 

t_flag=l; 

g=0; 

} 

else g+ +; 

/* Displays energy information gathered above. */ 

if(graphics_flag= = 1) display_graphicsO; 

if((text_flag= = l)&&(energy_flag= = 1)) { 

printf("\n Total Energy at time %lf is %lf',model_time,energy); 

energy_counter + + ; 

if(energy_counter= =20) { 

printf("\n\nStrike any key to continue..."); 

while(kbhitO = =0); 

energy_cou nter = 0 ; 

} 

} 

if(phase_flag!=0) { 

phasplotO; 

phase_counter+ + ; /* Used for Poincare sections */ 

} 

} 



2 . Upper Cutoff, Global Drive, Equation of Motion 
eqn motionO { 

/* The following loop calculates a round of velocities for one 
iteration of the model clock. This is where the equations 
of motion need to be located, if one wishes to change the 
model’s physics! */ 
int i,j,k; 

DOFOR(i,no_pendulums) 
if(coordinate[i] > 100) { 
retrieve_state0; 
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} 

DOFOR(i,no_pendulums) { 
if(i = =0){ 

acceleration = gamma[i] *(coordinate[i + 1 ] + coordinate! no_pendulums- 1 ) 
-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])+2*eta*cos(2*omega*model_time-delta))*coordinate[i] 
+ alpha*CUB(coordinate(i]); 

} 

else if(i= =(no_penduIums-l)) { 

accel erat ion = gamma! i] *(coord inate! i- 1 ] + coord i nate!0] 

-2*coordinate!i])-beta*momentum!i] 

-(SQR(omegaO!i])+2*eta*cos(2*omega*modeI_time-deIta))*coordinate!i] 
+ alpha*CUB(coordinate!i]); 

} 

else { 

acceleration = gamma!i] *(coordinate!i + 1 ] + coordinate!!- 1 ] 
-2*coordinate!i])-beta*momentum!i] 

-(SQR(omegaO!i])-t-2*eta*cos(2*omega*modeI_time-deIta))*coordinateIi] 

-l-alpha*CUB(coordinate!i]); 

} 

old_momentum!i] = momentum!i]; 

momentum!!] = momentum!!] -f- acceIerat!on*t!me_!nt; 

} /*DOFOR*/ 

DOFOR(i,no_penduIums) { 
oIdold_coordinate!i]=old_coord!nate!!]; 
oId_coord!nate!!]=coord!nate!!]; 

!f(p!nned_elements!i] = =0) { 

coordinate!!] = coordinate!!] -I- momentum!!] *time_int; 

} 



3. Lower Cutoff, End Driven, Equation of Motion 
eqn_mot!onO { 

/* The following loop calculates a round of velocities for one 
iteration of the model clock. This is where the equations 
of motion need to be located, if one wishes to change the 
model’s physics! */ 
int i,j,k; 

DOFOR(i , no_pendu 1 u ms) 
if(coordinate!i] > 1000) { 
retrieve_stateO; 

} 

DOFOR(i,no_pendulums) { 
if(i==0) { 

acceleration =gamma!0]*(coordinate!l]-coordinate!0]) 
-beta*momentum!0] 
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-((coordinate[i]/sqrt(l + 2*SQR(coordinate[i]))) +eta*cos(omega*model_time)); 

} 

else if(i= =no_penduIums-l){ 

acceleration =gamma[i]*(coordinate[i-l]-coordinate[i]) 

-beta*momentum[i] 

-(coordinate[i]/sqrt(l + 2*SQR(coordinate[i]))); 

} 

else { 

acceleration = gamma[i] *(coord inate[i + 1 ] + coordinate[i- 1 ] 
-2*coordinate[i])-beta*momentuni[i] 

*(coordinate[i]/sqrt(l + 2*SQR(coordinate[i]))); 

} 

old_monientum[i] = momentum[i] ; 

momentum[i] =momentum[i] +acceleration*time_int; 

} /* DOFOR */ 

DOFOR(i,no_pendulums) { 
old_coordinate[i] = coordinate[i] ; 
if(pinned_elements[i]= =0) { 
coordinate[i]=coordinate[i]+momentum[i]*time_int; 

} 



4. Upper Cutoff, End Driven, Equation of Motion 
eqn motionO { 

/* The following loop calculates a round of velocities for one 
iteration of the model clock. This is where the equations 
of motion need to be located, if one wishes to change the 
model’s physics! */ 
int i,j,k; 

DOFOR(i , no_pendulums) 
if(coordinate[i] > 100) { 
retrieve_stateO; 

} 

DOFOR(i,no_pendulums) { 
if(i = =0){ 

acceleration =gamma[0]*(coordinate[ l]-coordinate[0]) 

-beta*momentum[0] 

-(SQR(omega0[0])*coordinate[0] + eta*cos(omega*model_time)) + alpha*CUB(coordinate[Oj); 

} 

else if(i ==no_pendulums-l){ 
acceleration=gamma[i]*(coordinate[i-l]-coordinate[i]) 

-beta*momentum[i] 

-(SQR(omegaO[i])*coordinate[i])+alpha*CUB(coordinate[i]); 

} 

else { 

acceleration = gamma[i] *(coord inate[ i + 1 ] + coord inate[ i- 1 ] 
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-2*coordinate[i])-beta*momentum[i] 
-(SQR(omegaO[i])*coordinate[i]) + alpha*CUB(coordinate[i]); 

} 

oId_momentum[i] =momentum[i]; 

momentum[i] = momentum[i] + acceleration*time_int; 

} /* DOFOR */ 

DOFOR(i,no_penduIums) { 
old_coordinate[i] =coordinate[i]; 
if(pinned_elements[i] = =0) { 
coordinate[i] =coordinate[i] +momentum[i]*time_int; 

} 
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